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Introduction 


The Taylor experiment on Couette flow between coaxial circular cylinders 
has been the subject of numerous theoretical and experimental studies [1], 
This flow is rich in complex phenomena; so rich, in fact, that they are still 
being discovered [2], and our understanding of them is far from complete. In 
a typical Taylor experiment, the inner cylinder rotates with a constant 
angular velocity while the outer cylinder, along with the top and bottom 
walls, are kept at rest. The relevant geometric parameters are the radius 
ratio, which is the ratio of the radii of the inner and outer cylinders , and 
the aspect ratio, which is the ratio of the length of the annulus to its 
width. The dynamic parameter is the Reynolds number based on the angular 
velocity of the inner cylinder and the annulus width. The Taylor-Couette flow 
is strongly dependent on all of these parameters. The theory for the infinite 
aspect ratio case (which neglects end wall effects) and its correspondence to 
the experiments in cylinders of necessarily finite but large aspect ratio are 
reviewed in Di Prima and Swinney [1] and Di Prima [3]. Benjamin [4] has 
developed a rigorous qualitative theory for the existence of nonunique steady 
states for confined flows and their stability and transition with particular 
reference to finite-length Taylor-Couette flow. His predictions have been 
confirmed in a series of experiments [5] in cylinders of short length with 
fixed end plates. They have been further confirmed by the numerical results 
of Cliff e and Mullin [6] who discretized the steady Navier Stokes equations by 
a Galerkin finite element method and solved the resulting nonlinear algebraic 
equations by the pseudo arclength continuation method of Keller [7], 

Among other numerical studies of the finite-length Taylor-Couette flow 
problem, those of Alziary de Roquefort and Grillaud [8], and Neitzel [9] are 



- 2 - 


worth mentioning. Both investigations were based on the time dependent 
vorticity-stream function formulation along with the equation for the 
azimuthal velocity. They used a finite difference method to solve for the 
steady state by a time asymptotic approach. Their numerical results show the 
axial structure for small Reynolds number, the smooth development of the flow 
with rapid increase in vortex activity for Reynolds numbers in quasi-critical 
range, and multiple states for sufficiently large Reynolds number. The aspect 
ratio being large, their results are only pertinent to the exchange phenomena 
beyond the two-cell and four-cell interactions examined by Benjamin. 

No theoretical work on finite-length Taylor-Couette flow has incorporated 
the correct boundary conditions for fixed end walls. An idealized version of 
the end-wall boundary conditions is due to Schaeffer [10], He introduces a 
parameter alpha in the boundary conditions, with 0 < a < 1. The parameter 
a interpolates linearly between the two extreme cases: a = 0 corresponds to 

the infinite length problem which accommodates the Couette flow as an exact 
solution, while a = 1 corresponds to the finite length problem with the 
correct no-slip conditions being applied on the end walls. 

For o = 0, Blennerhasset and Hall [11] have considered the linear 
stability problem in the small gap limit, and the key result was that the two- 
cell primary flow changes into four-cell at the aspect ratio of approximately 
2.6. This should be compared with a value of roughly 3.7 observed by Benjamin 
for radius ratio of 0.615. Hall [12] has further derived the amplitude 
equations for this problem. An interesting feature of these amplitude 

equations is a quadratic term (absent in the infinite aspect ratio case) which 
can introduce hysteresis and soften bifurcations into smooth transitions. 


For small nonzero values of a , Schaeffer [10] used a Lyapunov- Schmidt 
reduction procedure and the methods of singularity theory to obtain results 
applicable to the exchange processes between 2m and 2m+2 cells, m > 2. Hall 
[12] has studied the two cell-four cell exchange problem using amplitude 
equations and a perturbation method for small values of alpha. 

The purpose of our continuing research effort is to solve the unsteady 
Navier-Stokes equations by a highly accurate spectral collocation method with 
a view to elucidate the underlying processes leading to laminar-turbulent 
transition in the Taylor-Couette flow. The present work is confined to the 
evolution of two-cell and single-cell Taylor-Couette flows with specific 
reference to the experiments of Benjamin and Mullin [13], Lucke, et al. [14] 
and Aitta, et al [15], The main result of these studies is a second order 
transition from a two-cell flow, symmetric under reflection about the mid- 
plane, to an asymmetric single-cell flow that ensues with increasing Reynolds 
number beyond a certain critical value. 


Splitting Scheme 

The incompressible Navier-Stokes equations for axisymmetric flow in a 
cylindrical geometry are, in conservation form: 
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and u, v, w are the r, 0 and y components of velocity, respectively. For the 
case where the outer cylinder and end plates are stationary, and the inner 
cylinder rotates, the boundary conditions are: 

u = v = w = 0 on the outer cylinder r = R Q , 

and at the top and bottom wall y = ± y^ 
u = w - 0 on the inner cylinder r = 
v = QR^, 

In the present calculations, the azimuthal velocity v is split into two 

parts, v = V. + v, where v, satisfies: 
b b 



= 0 


v, = 0 at r = R , and at y = ± y. 
b o l 


= ftIL at r = R. 

The quantity is computed and stored at the start of a calculation, 

which proceeds with the computation of u, v, and w at each time step. Note 
that these three velocity components satisfy homogeneous boundary conditions. 

A splitting method is employed to advance the solution from 
t n to t n+1 . Writing the Navier-Stokes equations in vector notation with u 
representing the velocity, (u, v, w) , we have: 
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u t + u»Vu = - 7P + vV ja in D, the annulus 


V* u = 0 


(5) 


ii = 0 on the boundary T 


n * 

the split scheme first advances u to an intermediate solution ii by 


solving: 


* * * 2 * 
u. t + u * Vu - vV _u 


( 6 ) 


* * 
u = g on I\ 


* * 

The intermediate boundary condition u_ = j[ will be discussed 

* n+1 

subsequently. Finally, the solution is advanced from u to u via: 


u f 1 = -7P n+1 


V-u n+1 = 0 


(7) 


A n+1 

n • u =0 on F 


where n is unit normal to the boundary r . Note that the final, "pressure 
correction" step by itself is an inviscid calculation, and is well-posed when 
only boundary conditions on the normal component of velocity are enforced. At 
the end of the full step there exists a non-zero tangential component of 
velocity on the boundary. The magnitude of this slip velocity can be reduced 
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by a proper choice of the intermediate boundary condition on u_ • Marcus [16] 

•nt 

has described the difficulties which arise from the use of ja = 0 as the 
intermediate boundary condition. The conditions used here are based on the 
work of Fortin, et al. [17]. Using backward Euler time discretization for 
Eq. 7 yields: 


n+1 * . . rjT) 

u = u - At VP 


( 8 ) 


and the slip velocity on the boundary is given by 


n+1 1 * ( * | . nt) n+l | % 

t • u Ip = T ’ V u | r - At VP | r ) 


(9) 


A A n+1 i / 

If x • u = 0 on the boundary, then t • u^ |p = 0(At) 

A 

where t is the unit tangent to the boundary T. However, if VP 


n+1 


is expanded in Taylor series about t = t : 


VP n+1 = VP n + At vp£ + 0(At 2 ) 


and the second term is approximated by 


At VP^ = (VP n - VP n X ) + 0(At 2 ) 


then Eq. 8 becomes 


x • u n+1 1 r = T * I p “At (2VP n - VP n h + O(At^). 



-7- 


3 

Thus the slip velocity may be reduced to 0(At ) through the intermediate 
boundary condition 


A * A n 1 
t • u | - t • At (2VP - VP n A ). 


Of course the boundary condition n • ii - 0 is retained. 

The pressure step is actually carried out in two parts. First, the 
divergence of the first of Eqs. 8 yields: 


„2 D n+l _ 1 „ * 

VP - -rr V*U 
At — 


where V • u 11 =0 is enforced. Then the velocities are updated using 


n+1 * A ^n+l 

u = u - At VP 


Note that this formulation requires a boundary condition for the 
pressure. This poses a problem, since there is no natural boundary condition 
for pressure. Deriving a condition by enforcing the normal momentum equation 
at the boundary is a questionable practice, as the equation need not hold on 
the boundary at the differential level of the equations. This inconsistency 
often produces explosive instability in a spectral code. Fortunately, the 
split scheme yields a self-consistent pressure condition, 


VP| r = o 


since both n • u and n • 


are zero on the boundary. The error involved 
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in this specification is, we believe, related to the overall splitting error 
of the scheme. 

Zang and Hussaini [18], have extensively investigated a related split 
scheme in which two coordinate directions were periodic and the third employed 
general boundary conditions. Comparison between split and unsplit codes using 
the same discretization yielded agreement to five decimal digits. However, 
they utilized a staggered grid for pressure in the non-periodic direction, 
obviating the need for a pressure boundary condition. Staggered grids in two 
dimensions either lack the ability to set both velocity components equal to 
zero on all boundaries, or are susceptible to an oscillatory "checkerboard" 
pressure mode. The unstaggered scheme used here, on the other hand, requires 
a pressure boundary condition. Actually, the consistent pressure equation 
derived by Zang and Hussaini yields exactly the same condition on pressure as 
used here. No instabilities were ever encountered with the present scheme. 

Discretization and Solution Scheme 

Since no-slip boundary conditions are enforced in both the r and y 
directions, Chebyshev spectral representation is appropriate in both 

directions. Collocation is used for a number of reasons: straightforward 

treatment of nonlinear terms and boundary conditions; capability to include 
coordinate stretchings; and ability to solve the resulting discrete equations 
rapidly. For further discussion of this form of discretization, see Hussaini 
et al. [11]. 

A coordinate stretching was employed in the radial direction to resolve 
the large gradients near the inner cylinder. The form of the stretching is: 
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[1 + b exp (-a)] (R - R , ) 

r = — - - p ^ r + R. (14) 

[1 + b exp (-a r^)] c l 

where r c is the radial coordinate in the computational space. Values of a = 2 
and b = 5 to 50 were typical in this work. The y-direction was not stretched. 

Time discretization of the first step in the split scheme involved the 
low-storage mixed Runge-Kutta/Crank-Nicholson scheme [20]. Writing the semi- 
discrete equation for the first step as 

u t = A(u) + D(u) (15) 

where A(u) and D(u) represent advection and diffusion terms, respectively, the 
mixed scheme advances from time step t n to t n+ * using a third-order Runge- 
Kutta scheme for the advection terms and Crank-Nicholson for the diffusion 
terms : 

u° = u(t n ) 

H 1 = At A(u° ) 


1 

u 


2 

u 


3 

u 


u° + J H 1 + i At (D(u° ) + D (u 1 )) 

H 2 = At A(ub - H 1 
ul + -fl H 2 + I 4 At (DCu 1 ) + D (u 2 )) 

H 3 = At A(u 2 ) - HI H 2 

u 2 + I 5 H 3 + i At (D(u 2 ) + D (u 3 )) 


( 16 ) 
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/.n+1 v 3 
u (t ) = u 


The second step of the split scheme uses backward Euler time 
discretization, as mentioned in the previous section. 

The above scheme is stable to 0(1) Courant numbers, which involves time 
steps many times larger than that desired for accuracy. Typically a time step 
which resulted in Courant numbers of .1 to .2 based on the smallest physical 
mesh spacing was used. The slip velocity resulting from the choice of time 
step was normally eight to ten orders of magnitude below the maximum velocity 
in a given direction. 

Note that the above scheme reduces the problem to a sequence of uncoupled 
Helmholtz/Poisson equations to advance the discrete solution. Since one time 
step requires the solution of nine positive-definite Helmholtz equations with 
Dirichlet boundary conditions, and one Poisson equation with pure Neumann 
boundary conditions, a computationally efficient technique had to be developed 
if this study is to be feasible. The present scheme is fairly efficient. 
This scheme involves preconditioned Richardson iteration, in which an optimum 
relaxation parameter is chosen dynamically using a principle of either 
minimum-residual (MR) or orthogonal-residual (OR) procedure. Details of this 
technique are as follows. 

Write the equation to be solved as 


L(u) = f 


(17) 


The residual at a given iterate "m" is 


R m = L(u m ) - f 


(18) 
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A preconditioned iteration scheme may be looked at in the following 
way. At a given iteration, the goal is to compute an update such that the 
next residual is zero, i.e. 


R m+1 = L(u m + Au) = 0 


(19) 


or for a linear operator L (as in this case): 


L(u m ) + L(Au) = R m + L(Au) = 0 


( 20 ) 


Approximating the operator L which is difficult to invert by a more easily- 
inverted operator M yields the following preconditioned scheme for the update: 

Au = M -1 [L(u m )] = M~ 1 [ R m ] (21) 


A well-known method for computing an optimal u) m involves minimizing the L 2 

j 

norm of the residual, R : 


• / t> 1 t* 1 \ t / _ m _ m v n tn / _ m _ / . r ^ j \ \ 

min (R , R ) = nan (R , R ) * 2w (R , L(Au)) 

ca 

2 

rn /v> 

+ to (L(Au) , L(Au) ) 

for which a stationary point is 

m _ (R™, L(Au) ) 

(0 = — 

(L(Au) , L(Au) ) 


( 22 ) 


( 23 ) 


This is the usual minimum residual (MR) method 
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An alternate method for chosing is to require that the successive 

residuals be orthogonal to each other. This yields an orthogonal residual (OR) 
scheme related to the method of steepest descent: 

(R m , R m+1 ) = (R m , R m ) + a) m (R m , L(Au) ) = 0 (24) 

and hence 

/ = - . ( . R — 2 — L J (25) 

(R®, L(Au) ) 

The scheme proceeds as follows from an initial guess u° and residual 
R° = L(u° ) - f: 


Au = M [R ] 
m+1 m . m ~ 

u = u + a) Au (2o) 

n m+l n m m _ * 

R = R - a) L(Au) 


etc. 

The preconditioning operator M is the second-order finite-difference 
operator which corresponds to the spectral operator L. Unequal-mesh spacing 
formulae, based on the physical-space point locations (which includes the 
radial stretching), are used in constructing M. A fixed V-cycle multigrid 
scheme driven by approximate-factorization is used to invert the 
preconditioning operator (Eq 21). It was found that complete convergence of 
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this step is not required for the overall preconditioned scheme to converge: 
essentially no difference in convergence rate is seen between fully solving 
the preconditioning step and merely reducing the residual (associated with 
inversion of the preconditioner) by two orders of magnitude which usually 
requires only two or three multigrid cycles. 

The overall convergence rate based on spectral operator evaluations is 
never slower than .35 (for Dirichlet boundary conditions); rates between .15 
and .2 are typical even on the highly stretched mesh. The use of MR or OR to 
compute the optimum relaxation factor yields equivalent convergence rates, 
although the OR scheme is found to be more robust in other contexts. 

However, the observed convergence rate deteriorates for this scheme when 
it is applied to the pure Neumann problem encountered in the Eqs (11) and 
(13). Since a large number of such problems are to be solved, a boundary 
influence-matrix technique is devised that reduces the computation time 
associated with this step. To compute the solution of 


L(u) = f in D 


u * 0 on T 
n 


a series of homogeneous solutions 



is obtained: 


L(u (i) ) = 0 in D 
u (i > 


(27) 


= on T 


( 28 ) 
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for each boundary point x^. For each solution , the vector of normal 
gradients [u^^ ] at all boundary points is computed. These vectors are 
collected as columns of a (square) matrix. This matrix is the influence of a 
unit disturbance at the boundary on the normal gradient at the boundary of the 
solution to the linear operator L(u) = 0. Denote this matrix by 

C = col {[u n (1) ], [u n (2) ], ...[u n (Nfe) ]} (29) 

where a total of points lie on the boundary. The matrix C may be inverted 
to yield the influence of a unit normal gradient at the boundary on the 
boundary value of a solution of L(u) = 0 , with provision for setting the 
level of the solution required. Then the solution of the original problem 
(Eq. 27) proceeds by solving the related problem with homogeneous Dirichlet 
boundary conditions: 


L(u ) = f in D 


0 on r 


(30) 


The vector of normal derivatives of u^ at the boundary [u^] is computed. 

u h must be corrected if it is to satisfy the desired Neumann boundary 
conditions. This correction is computed by applying the influence matrix to 
the boundary gradient error vector: 
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L(u C ) = 0 in D 

u C = C’ 1 [u h ] on r 
n 


(31) 


h c 

The desired solution is u - u + u • 

Computation and inversion of the influence matrix is done in a 
preprocessing step; since it is a function only of the mesh geometry, it may 
be stored and used whenever it is needed. To obtain a solution of the pure 
Neuraann-Poisson problem at each time step requires that we solve just two 
Dirichlet-Poisson problems and the solutions satisfy the desired boundary 
conditions exactly. 


Implementation and Performance : 

The above split scheme for solving the time-dependent incompressible 
Navier-Stokes equations has been inplemented with a view to processing on both 
scalar and vector computers. In scalar form on a CYBER-175, the code requires 
approximately 20 seconds per time step on a 17 x 17 mesh, and 56 seconds/step 
on a 25 x 25 mesh. A large fraction of these times is spent in computing the 
11 Helmholtz/Poisson solutions required in each time step. On the VPS-32 
vector processing machine at NASA Langley, the vector code requires about 14 
seconds step on a mesh comprised of 65 points in the y direction and 33 points 
in the r direction. Almost 85 percent of this time is taken in inverting the 
preconditioning step of the Helmholtz/Poisson solution as the approximate 
factorization part of the step is not vectorizable in a manner which yields 
adequate vector lengths for the two-dimensional problems. This observed 
performance for the two-dimensional Chebyshev-Chebyshev code is in line with 


the performance quoted in (21) for a three-dimensional Chebyshev-Fourier- 
Fourier code (both codes have similar operation counts). For a coordinate 
direction discretized with Chebyshev series using N points, the operation 
count for that direction is O(N^), whereas a Fourier-discretized direction 
requires only 0(N) operations. Thus, the present Chebyshev-Chebyshev method 
and the Chebyshev-Fourier-Fourier method both have total operation counts of 
O(N^), with the latter method having the advantage of larger possible vector 
lengths. 

When the azimuthal direction is added to this simulation, however, the 
relative performance improves dramatically. Since the azimuthal direction is 
periodic, Fourier series is an appropriate discretization. The computations 
are performed in the Fourier wave-space of that direction; the discrete 
equations for each azimuthal mode decouple (see [21] for details)), allowing 
vector lengths to increase by a factor of the number of points in the 
azimuthal direction. This increase in vector length improves the CPU 
seconds/point/ time step performance of the present scheme by a factor of four, 
and should allow the planned three-dimensional simulations to be performed. 


RESULTS 

The results presented here pertain to the axisymmetric two-cell/one-cell 
bifurcation, which occurs when the Taylor apparatus has an aspect ratio up to 
about 1.5. The form of the bifurcation depends sensitively on this parameter; 
experiments [15] show that this transition can change from supercritical to 
subcritical with variations in the aspect ratio of as little as eight percent. 


Most of the results are obtained either by making quasi-static changes in 
the Reynolds number and allowing the stable, dominant mode to settle, or by 
slowly sweeping through a Reynolds number range, monitoring the change in a 
particular mode* Of course, using a time-accurate code to simulate the 
bifurcations of steady-state solutions is quite inefficient, owing to the 
extremely small growth rates near the bifurcation points. However, the 
eventual aim of this work is to simulate the turbulence and broadband 
structure exhibited by Taylor-Couette flow at moderate Reynolds numbers and 
the code was developed with these time dependent flows in mind. The ability 
to simulate accurately the sensitive, steady state bifurcations at lower 
Reynolds numbers is an excellent test for the numerical method. 

Moreover, by using a time-dependent computation to investigate the steady 
state bifurcations, we can obtain information on the path which the system 
follows as states exchange stability. This is illustrated by the following 
result. A 17 x 17 grid is used in all these simulations with a few 
calculations on a 25 x 25 mesh for accuracy check. The time step in these 
simulations corresponds to a maximum Courant number of about 0.2; the time 
step is limited by accuracy, and not by stability. For the geometry of 
Benjamin and Mullin [13] with a radius ratio .615, and aspect ratio 1.05, the 
symmetric two-cell mode is allowed to stabilize at a relatively low Reynolds 
number (R = 62). The Reynolds number is then raised impulsively to 165, above 
the experimental bifurcation boundary of about 150 and the growth of the one- 
cell asymmetric mode is observed. Random machine roundoff error on the order 

-14 

of 10 provides the initial energy for the mode. The order parameter used 
here to quantify the asymmetry of the mode is due to Lucke et al* [14]: 
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. = [ drdz (u(r,z) - u(r,-z)) 

V /drdz ( | u(r, z) | + |u(r,-z)|) 

I 

The integrals were performed by spectral collocation. The logarithm of 

this parameter is shown in Fig. 1 plotted against time in units of the 

2 

diffusion time scale y /v . As can be seen, the initial instability leading 

Xf 

to the one-cell mode appears to be linear; that is, exponential growth with 
time is observed with only the later stages being modulated by nonlinear 
effects. Also shown in Fig. 1 is a plot of the logarithm of the disturbance 

i 

energy versus time. After an initial period, the disturbance energy grows at 1 

! 

a rate which is within 2% of double the growth rate of as expected. j 

Streamlines in a cross-sectional plane at various stages in the two- 
cell/one-cell exchange are shown in Fig. 2. Locations of these intermediate 
states on the ^ vs time curve are indicated in Fig. 1. Note that the 

progression between states is smooth without abrupt collapse or alteration of 
the flowfield structure. 

The geometry of Lucke, et al. [14] has a radius ratio of .5066, and an ■ 

I 

aspect ratio of 1.05. This choice of parameters leads to a smooth j 

! 

supercritical bifurcation to the one-cell mode as the Reynolds number is < 

increased beyond a critical value. A plot of the order parameter against 

Reynolds number for this simulation is shown in Fig. 3. This curve is traced 1 

in the direction of both increasing and decreasing Reynolds number. As we 

approach the initial critical value, the Reynolds number is varied very slowly 
at a rate of about ± .2 units based on the diffusion time scale. The growth 
rates are much larger on the upper branch of this plot, which permits larger 
changes in Reynolds number. This bifurcation diagram is validated by 
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restarting the simulation at various points along the curve and allowing the 
flow to settle to eight decimal digits for a fixed Reynolds number. 

Also shown in Fig, 3 are the results of Lucke et al. [14] for this 
case. Their results were computed using a staggered-grid finite-difference 
scheme on a grid of about 28 x 30 points; however, in that calculation the 
Reynolds number was changed at a rate of about 4,1 based on the present time 
scale, A large discrepancy in the critical Reynolds number as predicted by 
these two studies is noted. 

We also investigated the geometry of Aitta, et al, [15] where the radius 
ratio is .5, Their experimental results relate to three aspect ratios: 

1.129, 1.266, and 1.281; for which the two-cell/one-cell bifurcations are 

supercritical, transcritical, and subcritical, respectively. In Fig. 4, an 
order parameter due to Aitta et al. [15] is plotted as a function of Reynolds 
number. 

+y o 

J w(r,y) dy 



^ + y„ 

/ v | w (r ,y) | dy 

y i 

where r * R^ + .14(R 0 - R^) . Also shown in Fig. 4, are the experimental 
results of Aitta, et al. The two loci of states agree in form. They also 

agree as to the level at which the asymmetric branch becomes unstable and the 
width of the region in which both symmetric and asymmetric modes are stable. 
The critical Reynolds number from the simulation is within 3% of that of the 
experiment. The growth rates of the one-cell mode in the "hysteresis" region 
of the bifurcation are exceedingly small, several orders of magnitude below 
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those observed in the first two cases, producing a flow costly to 
accurately. Resolution requirements are also large; a 33 x 33 
required for this simulation. 


simulate 
mesh is 
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